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ABSTRACT 

We give an overview of state-of-the-art multi-loop Feynman 
diagram computations, and explain how we use symbolic ma- 
nipulation to generate renormalized integrals that are then 
evaluated numerically. We explain how we automate BPHZ 
renormalization using "henges" and "sectors", and give a 
brief description of the symbolic tensor and Dirac 7-matrix 
manipulation that is required. We shall compare the use 
of general computer algebra systems such as Maple with 
domain-specific languages such as FORM, highlighting in par- 
ticular memory management issues. 

Categories and Subject Descriptors 

J. 2 [Computer Applications]: Physical sciences and en- 
gineering — Physics; 1.1.2,4 [Computing Methodologies]: 

Symbolic and algebraic manipulation — Algorithms, Applica- 
tions 

General Terms 

Algorithms, Languages. 

Keywords 

Quantum Field Theory, Renormalization Theory, Feynman 
Diagrams. 

1. INTRODUCTION 

It has been saicQ that "in the thirties, under the demoral- 
izing influence of quantum-theoretic perturbation theory, the 
mathematics required of a theoretical physicist was reduced 
to a rudimentary knowledge of the Latin and Greek alpha- 
bets." Likewise, the ability to evaluate multi-loop Feynman 
diagram can be reduced to an ability to do some simple 
tensor algebra, some graphical manipulations to purge the 
calculations of divergences, a knowledge of some basic prop- 
erties of F functions, and the ability to evaluate integrals 
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numerically. Since this is obviously damaging to the delicate 
egos of theoretical physicists we advocate that the entire pro- 
cedure can — and should — be automated and handed over 
to computers which do not have egos and, as is well-known, 
never make mistakes. 

Precision measurements related to properties of elemen- 
tary particles are nowadays on a level which makes the in- 
clusion of quantum corrections to theoretical pre- and post- 
dictions mandatory. As will be explained below these quan- 
tum corrections are evaluated in the context of quantum field 
theoretic perturbation theory, the loop expansion. This al- 
lows for a systematic evaluation of quantum corrections to 
— in principle — any order in the coupling parameters of 
the theory under consideration. 

1.1 Loops and Legs: the State of the Art 

The complexity of perturbative calculations grows with 
the number of free internal momenta which are integrated 
over (loops), and the number of external particles of the pro- 
cess under consideration (legs). Only a limited number of 
observables have actually been calculated beyond the one- 
loop level. For one-loop scattering amplitudes on the other 
hand only one complete scattering process with more that 
five external legs has been evaluated up to now: the elec- 
troweak corrections to four fermion production in electron 
positron collisions which is highly relevant for e + e~-collider 
phenomenology 23 . For hadron colliders a lot of progress 
has been made very recently concerning the evaluation of 
six-point amplitudes. Different methods have been designed 
to tackle this problem and benchmark amplitudes for six- 
photon [STirTSirjH] and six-gluon scattering [T51 151 \M El [50] 
are now available. In the one-loop case so-called unitarity 
based methods have played a prominent role in these devel- 
opments (see [6] for a review and references therein). The 
unitarity approach is not based on Feynman diagrams but 
rather on bootstrapping tree level amplitudes. Although it 
has a big potential its applicability to amplitudes with inter- 
nal and external masses needs further developments before 
it may replace the Feynman diagrammatic approach. 

It should be noted that the closed-form structure of one- 
loop amplitudes is completely understood. One can show 
that each one-loop Feynman integral is expressible in terms 
of dilogarithms. Even at two loops the situation is less clear: 
only a limited number of four-point amplitudes, mainly for 
massless particles, are known. The corresponding Feynman 
diagrams are reduced first to an independent set of basis 
functions, the so-called Master integrals, by using integra- 
tion by parts identities [T7] and relations from Lorentz in- 
variance [26]. The analytic knowledge of the amplitude re- 



lies then on the successful evaluation of the Master integrals. 
This is a formidable task and typically only successful if the 
number of scales is very restricted. For example a break- 
through in that direction was the analytical evaluation based 
on Mellin-Barnes representations of the massless planar [411 
127] and non-planar [35] four-point two loop integrals in 2000. 
Although there has been some progress for more complicated 
cases it has been very slow. Unfortunately, even direct nu- 
merical evaluation is highly non-trivial as the presence of 
infrared (IR) and ultraviolet (UV) divergences within the 
graphs necessitates regularization; even in their absence in- 
tegrable singularities due to internal thresholds requires new 
numerical methods [71 l40l 151 l22l fT] . 

1.2 Beyond Two Loops 

Beyond the two loop level essentially only one-scale prob- 
lems have been evaluated so far. A very remarkable result 
here is the recent analytical evaluation of the three-loop 
splitting functions in QCD 49, 36] . By mapping the split- 
ting functions to Mellin moment integrals one can generate 
a hierarchy of difference equations which have to be solved 
recursively. The project took many person-years and needed 
to accumulate a database of integrals which occupied about 
3.5 GBytes of disk space. The algebraic manipulations relied 
on the domain-specific language FORM, or more precisely on 
constant improvements of this code. Indeed, the successful 
completion of this important calculation relied on the fact 
that the author of the algebraic manipulation code was a 
member of the collaboration. 

The FORM code was also the basis of the few known four- 
loop computations in perturbative QCD. All renormalization 
constants are known at this level [16]. An important mile- 
stone was in this respect the evaluation of the four-loop /3- 
function in QCD [46] in the so-called MS renormalization 
scheme. The problem was mapped to the evaluation of three 
loop propagator functions which were computed with the 
FORM package MINCER. The project amounted to the evalua- 
tion of about 50,000 Feynman diagrams and its confirmation 
by an independent group was an important issue 21 . The 
latter mapped the problem to the evaluation of the divergent 
part of some four-loop tadpole master integrals which were 
evaluated by solving integration by parts identities using a 
dedicated C++ code. To give a prominent example of a five- 
loop contribution to a very precisely measured observable 
we should mention the numerical evaluation of the anoma- 
lous magnetic moment of leptons by Kinoshita et. al. Some 
dominant contributions are known numerically now even at 
the five-loop level [33] . The parameter integrals correspond- 
ing to the Feynman diagrams are evaluated with the Monte 
Carlo integration routine Vegas. This heroic effort will pre- 
sumably not be matched by an analytic computation in the 
foreseeable future. Remarkably enough the three loop fully 
analytic result is available 35 . 



2. QUANTUM FIELD THEORY 

When quantum field theories (QFTs) were developed in 
the 1930s they were expressed in terms of field operators 
acting on Hilbert spaces and satisfying various commutation 
and anticommutation relations. Nowadays we prefer to think 
of them in terms of functional (or path) integrals, which give 
us greater geometrical insight (for example for the introduc- 
tion of Fade'ev-Popov ghosts for gauge-fixing non-abelian 



gauge theories) with just as high a level of mathematical 
rigour (namely not very much). 

2.1 QFT as a Functional Integral 

The quantities we want to calculate are matrix elements 
of the scattering (S) matrix, and these can be written as 
integrals over all possible field configurations with a measure 
exp(iS ((/>)) where S = J d n xL((j),d(j)) is the action which is 
the space-time integral of the Lagrangian density L defining 
the theory. We take the dimension of space-time to be n 
rather than four as this is a convenient way of "regulating" 
the divergences that occur; when we have removed all the 
divergences by some "renormalization" procedure that we 
shall discuss later then we will take the limit that n — > 4 
(or n — > 6 which is convenient for the scalar (f> 3 theory we 
are using as an example). For a scalar field theory, which 
for simplicity is all we shall consider her^3 the Lagrangian 
is L = \4>{d 2 + m 2 )4> + A<^) 3 . A typical S'-matrix element is 
then 



{t/>(x)4(ji)\S\4>(s)) = ± J d<fi<j>(x)<fi(y)<l>(z 



iS(<t>) 



where the points x and y are at time t = oo and z is at t = 
-co. The integral is over all fields <f>: that is d<f> = Yl x d<f>(x) 
where the product is taken over all space-time points. We 
shall not dwell on the precise definition of such integrals but 
we will consider the situation where the coupling constant 
A is small and we can apply perturbation theory. Formally 
one can consider that the weight factor exp(iS') should be 
written as exp(iS/h) in order to get the dimensions correct, 
and that as Planck's constant h is very small the integral 
will be dominated by field values near to the minimum of the 
action which occurs at <j> — 0. In reality this is completely 
bogus, firstly because we are interested in phenomena where 
the natural scale is set by working in units where h = c = 1, 
and secondly because our cj) 3 action is not bounded below 
as (f> — > — oo. To justify our argument we need to point 
out that we are really expanding in the powers of A rather 
than h, and that for small perturbations about <f> = the 
theory does not notice the "non-perturbative" instability of 
the vacuum. We could try to be a little less cavalier and 
consider a 4 interaction instead of or as well as the <f) 3 one, 
but it turns out that this theory is also sick, and so we will 
keep things simple and just ignore these problems. 

2.2 Perturbation Theory 

The perturbative expansion we use is just an (unjustified) 
generalization of the method of steepest descents to an infi- 
nite dimensional integral. We add a linear source term J to 
the action so that we can write 



/ 



d(f> <j){x)<j)(y)<t){z)e 



iS(4>) 



H) 3 



d 3 Z(J; A) 



dJ{x)dJ{y)dJ{z) 



where 



Z(J; A) 



exp 



i J d n x , 



dJ{x) 



Z (J; 0) 



where we have written J<j) for f d n x J(x)cj>(x), and Zo(J) = 
/#e i (i w+,/ '*) oc e-i JAJ , where we have introduced the 

2 The reader should trust us that realistic theories with gauge 
fields, fermions, ghosts and so forth just add some technical 
complications but do not essentially alter the strategy and 
methods we propose. 



kernel K = d 2 + m 2 and its Green's function A = K~ Y 
with appropriate boundary conditions. We obtain the re- 
quired asymptotic expansion by expanding exp (iA9 3 /<9J 3 ) 
as a Taylor series, and then each resulting term can be drawn 
as a Feynman diagram with the rules that there is a factor of 
iX associated with each vertex whose position is integrated 
over all space-time locations, and a propagator iA with each 
edge. If we Fourier transform to momentum space then we 
obtain the equivalent diagrams but with the rules that mo- 
mentum is conserved at each vertex and there is an integral 
over the momentum flowing round each loop, and the prop- 
agators become A(fc) = (k 2 — m 2 + i0+) _1 . 

3. DIVERGENCES AND 
RENORMALIZATION 

The infamous UV divergence disease of QFT is now im- 
mediately apparent, it manifests itself even if we work in 
Euclidean rather than Minkowski space. Even the simplest 
one-loop diagram corresponding to the Euclidean integral 

r <r_k r dfciifcH"- 1 

[P '~J [k 2 +m 2 ][(k + p) 2 +m 2 ] <C J max(||fc||,m) 4 

for some constant c, and this clearly diverges in n > 4 di- 
mensions. 

What saves the day is that the divergence is local (a Dirac 
5-function or its derivatives in space-time, or equivalently 
a polynomial in the external momentum p in momentum 
space). We may see this easily by differentiating the Feyn- 
man integral with respect to p: 

dl(p) = _ f d n k2(k + p) fl 

dpv ~ J [k 2 + m?][{k+p) 2 + m 2 ] 2 ' 

For large ||fc|| the integrand now behaves as ||fc|| -5 rather 
than ||fe|| -4 as before; differentiating d = n — 3 times suffices 
to render the integral convergent. We may therefore express 
the Feynman integral as 

up) = T d i( P ) + r d n P1 . . ■ r^ 1 drp d d d i( Pd ) 

Jo Jo 

by iterating the fundamental theorem of calculus d times 
(otherwise known as Taylor's theorem) where T d I(p) is a 
polynomial of degree d with divergent coefficients. We have 
simplified the notation by writing d for 9/<9p M , suppressing 
all the tensor indices (such as fi), and lumping together all 
the external momenta into one p. 

We can remove the offending divergent polynomial by add- 
ing it to the action as a new term giving rise to a vertex 
that is formally of one-loop order. The number of loops in 
a Feynman diagram is formally equal to the corresponding 
power of h in the perturbation expansion, so we may choose 
to imagine that all the couplings in the action, A, m, and 
the coefficient Z of the kinetic term (f>d 2 <f> that we forgot to 
include originally, may be expanded in powers of h where all 
the terms but the first just serve to remove the unwanted 
divergences. Clearly there are two requirements that must 
be satisfied for this to work: (i) all the divergences must be 
of the form of terms that occur in the action, and (ii) all 
the divergences must add up in just the right way to cor- 
respond to such counterterms. We usually rephrase (i) by 
saying that we must include all monomials in the action that 
are allowed (i.e., which do not violate physically necessary 
conditions such as locality and unitarity), as if we do not 



we will have to come back later and put them in so as to 
cancel the divergences that arise. Of course, we can only 
allow there to be a finite number of parameters (and thus 
monomials) in the action if we are to have a theory with any 
predictive power, but happily we can get away with includ- 
ing only monomials of dimension less than n, as this suffices 
to cancel all possible divergences by a simple power-counting 
argument (q.v., H3.5p . 

3.1 Regulators 

Of course this is all rather messy and embarrassing, so we 
reformulate the procedure in the following language: in order 
to define a QFT we need to introduce a regulator of some 
kind which makes our manipulations mathematically well- 
defined. We then renormalize the theory by adjusting the 
"bare" coefficients in the action to absorb all the would-be 
divergences, and finally we take the regulator away to ob- 
tain a well-defined finite theory as the limit of a renormalized 
regularized functional integral. There are many choices of 
regulator, working in n dimensions ("dimensional regular- 
ization" 44, 431 1121 [18] ) is just one rather convenient pos- 
sibility, and they are all supposed to give the same answer. 
What turns out to be crucial is that the regulator preserves 
as many of the symmetries of the original theory as possible, 
as we can then exclude counterterms that do not also have 
these symmetries. 

3.2 Renormalization Conditions and 
the Renormalization Group 

In order to fix the finite parts of the counterterms in the 
action we need to specify a set of renormalization conditions, 
i.e., a set of experimentally determined values for some pro- 
cesses that can be solved for the parameters in the action. 
This is no different from what happens in classical theories, 
except that the parameters themselves have a less obvious 
physical meaning. Indeed, we can choose the renormaliz- 
ation conditions in many different ways; for example we can 
specify them as a function of some energy scale fj,, and the 
fact that they all fix the same physical theory tells us the 
physical quantities must not depend upon p. The invari- 
ance under this group of reparameterizations is known as 
renormalization group invariance, and is very useful because 
it is not in general respected by the perturbative expansion 
to any given order. 

The parameters in the action that are just constants on 
the classical level thus become in general scale dependent if 
radiative corrections are included. In the context of Quan- 
tum Chromodynamics (QCD), our theory of the strong in- 
teractions, this leads to asymptotic freedom, the logarithmic 
decrease of the coupling strength with increasing interaction 
energy. For QCD the perturbative approach is limited to the 
high energy regime where the gauge coupling parameter, a a 
is sufficiently small. 

In the case of gauge theories (which are the theories we 
really use to describe nature) we must choose renormaliz- 
ation conditions that preserve the gauge symmetries or the 
whole theory falls apart. On the other hand there are cer- 
tain symmetries, such as some forms of chiral symmetry (a 
peculiar symmetry that occurs because our theories involve 
fields that transform as spinor representations of the cover- 
ing group Spin(3, 1) of the Lorentz group SO(3, 1)) and scale 
invariance symmetry cannot be maintained by any regula- 
tor and thus are not symmetries of the quantum theory even 



though they look like perfectly good symmetries of the un- 
derlying classical action. There are physical implications of 
these so-called anomalies, such as the decay of a pion into 
two photons 7r — > 77, and this is evidence that renormaliz- 
ation is necessary and not just an ugly contrivance that we 
could avoid by being more clever. 

What will concern us for the rest of this paper is a more 
detailed investigation of (ii), namely showing that all the 
divergences are local and appear in the right way to be can- 
celled by counterterms in the action. For instance, if we go 
beyond one-loop Feynman diagrams are all the divergences 
still local? The trouble is that multi-loop Feynman diagrams 
are very complicated integrals, and it is hard to see what is 
going on without getting lost. Nevertheless, with a little ef- 



fort one can see that a typical two loop diagram 



has 



non-local divergences; however, before abandoning all hope 
we notice that the one-loop counterterms we introduced to 
cancel the one-loop divergences have to appear in loop dia- 
grams themselves. Their divergent coefficients thus multiply 
the non-local finite parts of the corresponding Feynman in- 
tegrals also leading to non-local divergences. What we must 
show is that these one loop countergraphs, which are for- 
mally of two loop order (0(ft 2 )) exactly cancel the non-local 
part of our two loop graphs, leaving only a local divergence 
that can be absorbed by adding local two loop counterterms 
to the action. 

3.3 Henges and the R and R Operations 

This is most easily done in two stages. The first is we de- 
fine a procedure that removes all divergences by local Taylor 
series subtractions, and then we show that the subtractions 
thus made add up to give local counterterms. We will discuss 
the first stage here, as it is what needs to be implemented to 
automate the renormalization process; the second is a purely 
combinatorial proof which the interested reader can find in 
the literature [131121119]. 



The two-loop diagram 



illustrates the difficulty we 



face. Clearly when all the loop momenta get large simulta- 



neously the diagram diverges as \\k\\ 



, and this overall 



divergence is local. However, what happens when one of 
the two three-line loop momenta gets large while the other 
stays small? How do we disentangling these overlapping sub- 
divergences? Our approach is a systematic decomposition of 
the space of all loop momenta based on the structure of the 
graph itself. 

We need only consider graphs which are one-particle ir- 
reducible (1PI), that is ones which remain connected when 
any line is cut. For any line £ in a, 1PI Feynman diagram Q 
we can decompose the graph uniquely into a single loop con- 
taining £ stringing together a set of 1PI subgraphs that we 
will call a Henge TL{TL,£) |15l 132] , We will discuss efficient 
ways of representing and computing Henges in Sj5] 

Clearly at any point the space of all loop momenta some 
line must carry the smallest (Euclidean) momentum^ Hence 
we may decompose an arbitrary Feynman diagram Q into a 
sum of contributions each from a region where a particular 
line I £ G carries the smallest momentum, and in each such 
region we can use our Henge decomposition to write the 



Feynman integral as 

h(G) = ^2j°° d n k e i ke (Q/n(g,£)) n h e (0), 



een(g,e 



where is the Feynman integral for the graph O re- 

stricted to the region of momentum space where all the lines 
carry momenta of magnitude greater than A, and i\(G/TC) is 
the product of all the lines in the single loop Q /Ti obtained 
by shrinking all the graphs in TL to points, again restricted 
to have momentum of magnitude greater than A. The "small 
momentum cutoff" A serves as a convenient parameter for 
our recursive definition. 

As an example let us consider the overlapping two-loop 
diagram we considered before. Its concomitant Feynman 
integral is 

I(p) = J d n kd n k' A(fe)A(fc + p)A(fe - fe')A(fe' +p)A(k'). 

The Henges that arise are i ' ' w ^ ere ^ e 

heavy lines show the Henges corresponding to each of the 
light lines. For example, the contribution from the region 
where the light line in the first diagram is smallest is 




d n k" A x (k")I k „ 



<D 



where A\(k) = A(fc)#(||fc|| — A) and the inner integral (cor- 
responding to the blob in the diagram on the left and the 
heavy loop in that on the right) is 

d n k A k „ (k)A k „ (k + p)A k „ (k + k" + p)A k „ (k + k"). 

life" II 

Observe how the graphical structure dictated the change of 
variable k" — k — k! to the loop momentum of the shrunken 
graph. 

With this decomposition it is obvious that the divergences 
can be put into two classes, overall divergences that occur 
when all the momenta simultaneously get large, and sub- 
divergences that are isolated to the subgraphs occurring in 
the Henges. Let us introduce an operator R that removes 
all the divergences from a given Feynman diagram: it first 
removes all the subdivergences (an operation that is called 
R) by recursively applying R to the elements of the Henges 



Rh(G) = J2 

ee5 



d n k l i k Ag/H{g,i)) Yl Mk t (P) 



een{g,i 



Up to a set of measure zero. 



and then it removes the overall divergence (if necessary) by 
subtracting the leading terms of the Taylor series expansion 
in the external momenta as described before for the one-loop 
case, RI X (G) = Ix(G)~T dcs{e) RI (G)- The number of terms 
removed by the Taylor subtraction operation T is fixed by 
simple power counting rules (q.v., 5j33}. Note the subtle but 
vital fact that the subtraction term has its small momentum 
cutoff set to zero. 

3.4 BPH and Zimmermann's Forests 

In order for this to be a valid renormalization procedure we 
need to prove two things: first that all the overall divergences 
are local — that is polynomial in the external momenta — 



so that they get eaten by the Taylor series subtraction, and 
second that all the subtractions that are made add up in such 
a way as to correspond to counterterms in the action. The 
former condition is equivalent to showing that the deriva- 
tives d des ^ +1 RI \(Q) with respect to the external momenta 
are finite, and follows by induction from (a) the inductive as- 
sumption that |7a(S)| < cmax(A, m) des< - e '' for some constant 
c, (b) the requirement that differentiation lowers the power- 
counting degree, deg(dQ) < deg(tj) — 1, and (c) the fact that 
differentiation commutes with the R operation, [d, R] — 0. 
To prove (c), as well as the condition (ii) above, we rewrite 
the R operation in purely graphical form (as it was originally 
defined by Bogoliubov and Parasiulfl), namely 

ri(q) = HG/s) * n (-r dcs(r) i?/(r)) . 

sew(S) res 

Here S is a spinney, that is a covering of Q by a set of 
1PI subgraphs (possibly including single vertices), and the 
proper wood VV is the set of all such spinneys (excluding the 
one consisting of just Q itself). The notation I(G/F) * /(F) 
means that the 1PI subgraph V in Q is to be shrunk to a 
point and replaced with the value /(F). 

Proving the equivalence of this definition of ii to the pre- 
vious one for A = is left as an exercise for the dedicated 
reader who need just show that in the recursive Henge de- 
composition the sequence of subtractions made always cor- 
respond to some spinney, and that for each spinney S the 
corresponding subtractions occur in the Henge decomposi- 
tion for every ordering of the line momenta in Q/S. 

One can expand out the recursion even more and obtain 
an explicit formula for all the subtractions that are made in 
a graph in terms of nested non-overlapping Taylor subtract- 
ions; the necessary graphical apparatus for this was given 
by Zimmermann |51j in his forest formula. Zimmermann 
introduced his forest formula because it explicitly defines 
a convergent integral corresponding to a Feynman diagram 
without the need for any intermediate regularization scheme. 
This property is the reason for our practical use of the R op- 
eration: we use it to generate the integrand of a manifestly 
convergent integral, although in practice this is more easily 
implemented using Henge recursion. From a formal point 
of view the utility of Zimmermann's approach is limited by 
the fact that although it defines finite integrals these are 
not directly related to the original functional integral in a 
straightforward way, and therefore the various identities be- 
tween Feynman integrals that follow because of this that 
are needed to guarantee desirable properties such as gauge 
invariance (Ward or BRS identities) or unitarity (cutting re- 
lations) are no longer manifest. Indeed the former are still 
untrue in the case of anomalies despite the fact that no reg- 
ulator is introduced. 

The fact that Bogoliubov's definition is purely graphical 
makes it obvious that [d, R] = 0, and it is also possible to 
give a combinatorial argument that shows that all the sub- 
tractions made, summed over all Feynman graphs, formally 
corresponds to adding —TRe 1 ^ to the action [T4l [2], 

3.5 Power Counting 

The BPH theorem tells us how to remove all the diver- 
gences from a Feynman diagram by local subtractions, but 

4 Their original proof was corr ected later by Hepp, hence it 
is known as the BPH theorem [ill E3 [28] [3] 1141132] 



it does not guarantee that the number of counterterms re- 
quired is finite. If it is not then the resulting theory has little 
if any predictive power, as an infinite number of renormal- 
ization conditions will be needed to fix the finite parts of all 
these new vertices in the theory. The class of theories that 
can be renormalized with a finite number of counterterms, 
called renormalizable theories, is therefore of central inter- 
est. We can give a simple power counting rule that tells us 
when a theory is renormalizable (but not the converse). 

The deg function introduced in i)3.3l must be chosen such 
that [A(fc)| < c ■ max(fc, m) deg ^ A ' for all vertices, propaga- 
tors (lines), and their derivatives in the Feynman rules as- 
sociating integrals with graphs; this is the starting point for 
the inductive proof that |/a((7)| < c ■ max(A, rn) deB ^ s ' that 
is used to establish the BPH theorem in our approach. For 
our working example of (f> 3 theory in n dimensions the degree 
of each propagator is d — —2 and the vertices have degree 
d! = 0. The degree of a graph with I internal lines, V ver- 
tices, and L loops is thus deg = Id + Vd' + Ln, recalling that 
there is an n-dimensional momentum integration associated 
with each loop. Now, any connected graph with V vertices 
has a spanning tree containing exactly V — 1 edges, and the 
remaining I — V + 1 edges each give rise to an independent 
loop; furthermore each internal line has two ends and, in our 
theory, each vertex connects three lines, so "conservation of 
line ends" implies that 3V = 21 + E where E is the number 
of external lines. Eliminating L and I from these equations 
we find that deg = n - E\(d + n) + V|(2d' + 3d + n) = 
n — Edim((/>) + Vdim(V) where the dimension of the field 
dim(</>) = |(n — 2) and of the vertex is dim^t^) = |(n — 6). 

These dimensions can be read directly from the Lagrangian 
describing the theory by some simple rules. We see that if 
dim(V) > then we can find arbitrarily overall divergent 
graphs just by increasing the number of vertices sufficiently, 
whereas if dim(V) < then only a finite number of graphs 
can have overall divergences (such theories are called super- 
renormalizable) . The case dim(V) = (renormalizable) 
is particularly interesting as there can be an infinitude of 
overall divergent graphs, but as their degree is bounded by 
Edim((j)) (provided this is positive) only a finite number of 
Green's functions are overall divergent. In four dimensions 
n = 4 we have dim(0) = 1 and dim(V) = — 1 for <f> 3 the- 
ory, so it is superrenorrnalizable, whereas in six dimensions 
dim(0) = 2 and dim(V) = it is renormalizable. 

4. NUMERICAL VERSUS ANALYTIC 
EVALUATION OF INTEGRALS 

Heretofore people have usually calculated Feynman inte- 
grals by using some regulator to make all the manipulations 
well-defined and then arranged the resulting expressions to 
cancel the would-be divergences between graphs and coun- 
tergraphs, leaving answers with a finite limit as the regula- 
tor is removed. This procedure is effective if all the integrals 
can be evaluated in closed form0 but it is very difficult to 
completely automate the procedure as the closed-form eval- 
uation of the integrals is a new challenge at each order of 
the loop expansion. 

5 In fact, what is really needed is that the would-be divergent 
parts can be evaluated in closed form. The usual folklore 
is that the work required to evaluate the divergent parts 
of diagrams with £ + 1 loops is about the same as that to 
evaluate the finite parts with £ loops. 



We wish to use the R operation to generate manifestly 
finite Feynman integrals, which we can then evaluate nu- 
merically. The advantages of this approach are obvious, but 
there is a price to pay. First, we need to evaluate the albeit 
finite integrals for every value of the external momenta and 
parameters (renormalization conditions) separately. Second, 
the integrals are over a unit cube in I — 1 dimensions (q.v., 
EJOJ where I is the number of edges in the graph, and thus we 
have to resort to Monte Carlo methods in all but the most 
trivial cases, so we are only able to achieve limited precision. 
Third, in most cases where perturbation theory is applicable 
we are expanding in a small parameter (the fine structure 
constant a « in QED for example), so there is no point 
in computing the £ + 1 loop contributions unless the errors 
at £ loops are sufficiently small. Our window of opportunity 
is to compute one loop higher order that the current best 
closed-form solutions. 



5. GRAPHICAL ALGORITHMS 

We now consider how we can implement the R operation 
in practice in our symbolic-numeric scheme. It is much eas- 
ier to implement the more recursive Henge approach than 
Bogoliubov's definition or Zimmermann's explicit forest for- 
mula, as the complexity of the problem is then naturally 
handled by a correspondingly recursive program. To do this 
we need some reasonably efficient graph manipulation al- 
gorithms, as the text-book methods using adjacency ma- 
trices and suchlike are slow in practice. Fortunately such 
algorithms are well-known, such as the Galler-Fisher [251 
134] algorithm for finding equivalence classes given pairwise 
equivalence relations, which may be equivalently viewed as 
a means of splitting a graph into its connected components. 

The basic idea is to represent a graph as a set of edges 
each represented as a pair of vertices, and to label the ver- 



tices with small integers; for instance the graph 



could 



be represented as {(12), (23), (34), (41), (24)}. We create an 
array F indexed by vertices and initialised to some special 
value (such as zero), and then run once through the edges 
(ij) computing the ancestor of each of its vertices, where the 
ancestor A(i) is % if F[i] = or otherwise. If the an- 

cestors are unequal A(i) ^ A(j) then we set F[A(i)] <— A(j). 
When we have exhausted the supply of edges the ancestor of 
any vertex labels the connected component that the vertex 
belongs to. 

This algorithm is pleasantly easy to extend to carry out 
other graphical tasks efficiently. To find all the irreducible 
edges in a connected graph (those that may be cut with- 
out disconnecting it) we just keep track of the route from 
a vertex to its ancestor by introducing an Z3 chairQ- valued 
array P indexed by vertices that contains the chain of edges 
connecting i to F[i], and defining the route R(i) from i to 
A(i) to be zero if i = A{i) or the chain R(i) = P[i]®R(F[i]) 
otherwise. Upon examining the edge (ij) we note that the 
chain Cy = R(j) 95 (ij) Q R(i) connects A(i) to A(j), so if 
A(i) ^ A(j) as well as updating F[i] as above we also set 
P[i] <— Cij. If, on the other hand, A(i) = A(j) then we have 
discovered that Cij is a closed loop, and all the edges occur- 
ring in this loop are marked as irreducible. All the edges not 
so marked by the end of the loop over edges are reducible. 



One use of this algorithm is to route momenta through 
a graph. We multiply each loop by a symbol for the cor- 
responding loop momentum, and for each vertex i at which 
external momentum p enters the graph we multiply the route 
R(i) by p. The coefficient of any edge in the sum of all these 
p-chainqj is the momentum flowing through that edge. 

Indeed, it is but a small step further to devise an efficient 
Henge-finding algorithm. Start with a 1PI graph Q and re- 
move any edge £; as Q is just a set of edges this graph is just 
the set of edges Q — {£}. Now apply the previous algorithm 
to this graph, partitioning the lines into two sets Q = I U R 
of irreducible lines and reducible lines, and apply the first 
algorithm to partition / into disconnected pieces / = {#}. 
The required Henge Ti(Q,£) is just this set, and the single 
loop g/H{g,£) =RU {£}. 

Using these algorithms we have implemented the subtract- 
ions specified by the R operation on an arbitrary diagram 
in 3 theory, and the generalization to more interesting the- 
ories with more complicated actions is both straightforward 
and in progress. It is perhaps worth noting for the benefit of 
the dedicated reader who has thought about the equivalence 
of the Henge and Bogoliubov definitions of R that too naive 
an implementation of the former will produce the same spin- 
ney multiple times, corresponding to different orderings of 
lines in Q — S. 



6. FEYNMAN PARAMETERS 

We have described in some detail the means of obtaining 
a finite Feynman integral, but we still have to evaluate it. 
There are various ways of doing this: we could for exam- 
ple directly evaluate the loop integrals in momentum space, 
leading to an nL- fold real integral. For some regulators, such 
as the lattice, where we discretize the functional integral by 
approximating space-time by a hypercubic grid, this is es- 
sentially the only way to proceed. In practice, however, we 
usual carry out perturbative QFT calculations using dimen- 
sional regularization in which case all of the propagators are 
essentially inverse quadratic forms in the momenta, as we 
used above. In this case there is a convenient F function 
identity that allows us to write a product of such propaga- 
tors as a power of a single quadratic form 



F(a)F(/3) 



= T(a + /3) 



dtt - 1 



A^B? v Jo (A + Bt) a +P' 

This can be iterated to establish the formula 



n 



(N-l)l / dx 



5 1- 



(i) 



the quantities x\ being known as Feynman parameters. 

Introducing Feynman parameters and interchanging the 
order of the momentum and parameter integrals (which is 
valid in the presence of a regulator) we can combine all nL — 
2ui loop momenta into a single vector and carry out the 
momentum integration using another V function identity, 



d^k 



k 2 +F(p,m) 



7r " r r( a r ) f(p ' m) " a 



A Z3 chain is a sum of edges with coefficients ±1 or 0. 



7 A p-chain consists of sums of edges with coefficients being 
± symbolic names for momenta. 



and related formulae obtained from this by differenting with 
respect to the external momentum p. We are then left with 
/ parameter integrals to evaluate, in our case numerically. 

6.1 Sectors 

It is interesting to see how the Henge decomposition is re- 
flected in Feynman parameter space. To this end we define 
a sector of Feynman parameter space as a region where the 
parameters are totally ordered, for example one such sec- 
tor is x\ > X2 > ■ ■ ■ > xn- If we restrict the integral in 
equation ([!]) to this sector we find 



dxi 



d,X2 ■ ■ 



r x N-i 
/ dxf 
Jo 



6 1- 



J2k=l x k) 



[E^Li x iQi] 

1 N 1 



so for example the contribution to the product Tuy rep- 
resented using Feynman parameters t, it, and v from the 



sector t > u > v is 



This means that the 



T(T+U)(T+U+V) ■ 

quadratic form denominator T corresponding to parameter 
t is guaranteed to be smaller than all the other denomina- 
tors, which means that it corresponds to a line carrying the 
smallest momentum in the Henge decomposition. 

We note in passing that especially in the context of over- 
lapping IR divergences the the iterated application of sector 
decomposition has proved to be a powerful way to deal nu- 
merically not only with multi-loop Feynman diagrams [7] [8] 
but also high-dimensional phase-space integrals [9]. 

7. TENSOR AND 7-MATRIX 
MANIPULATION 

In this pedagogical presentation we have illustrated our 
calculations using only scalar fields; in general we need to 
introduce vector fields for gauge bosons and spinor fields 
for fermions, these lead to the additional complication of 
having tensor expression in the Feynman integrand numer- 
ators. Actually our symbolic manipulation code already has 
to deal with such numerators, as they occur when we Tay- 
lor subtract even for scalar theories. What we end up with 
is some complicated SO(4) tensor expressions which need 
to be simplified into some canonical form; these expressions 
can themselves be viewed as graphs, where the tensors are 
the vertices and the edges are pairs of contracted indices, 
and we use the algorithms of section [J5] to carry out this 
simplification. 

For spinors, which strictly speaking are representations of 
the universal covering group Spin(n) of SO(n), the numer- 
ators also become functions of spinor operators built out of 
the Dirac (Clifford) algebra generators 7 M that satisfy the 
anticommutation relations {7 M ,7„} = 2g^ v where g is the 
metric tensor. Efficient algorithms are known for simpli- 
fying 7 matrix expressions 30, 31, 20 and have been im- 
plemented in the REDUCE package CVIT [29]; these algo- 
rithms are based on reducing the "spin-network" or "bird- 
track" diagrams in terms of sums of irreducible represent- 
ations of Spin(n) and their 3j and 6j coefficients. We are also 
currently investigating the generalization of these methods 
to handle arbitrary representations of SO(n), which might 
lead to much more efficient algorithms for symbolic tensor 
manipulation than exist at present. 



8. INFRARED CANCELLATIONS 

The dedicated reader might have observed that we do not 
live in Euclidean space, but in Minkowski space0 Fortu- 
nately this has absolutely no effect on the algebraic manip- 
ulations we have described: we end up with a renormalized 
integrand depending on a set of SO(4, C) invariants con- 
structed out of the external momenta, such as p 2 , p ■ q, 
£ l ±vpcrP' 1 q v r p s' T and the difference between the real forms 
SO(4, R) and SO(3, 1,R) is just that p 2 can be negative, 
for example. 

While the cancellation of UV divergences works just as 
well in Minkowski space as it does in Euclidean space, new 
IR divergences can occur. Previously we evaded all IR di- 
vergences by keeping the mass m of our particles non-zero, 
in the real world we need to handle the fact that some real 
particles are massless (such as the photorfl) and also that 
there are a variety of other IR singularities possible because 
k 2 = 7^ k = in Minkowski space. 

There are two different situations to consider: 

• There are true IR divergences that make the theory 
meaningless. For example, massless scalar field theory 
does not exist because it is IR sick in this manner. 

• There are IR singularities that occur at special val- 
ues of the external momenta. The emission of a zero- 
energy photon in QED exemplifies this situation. 

We shall reject theories in the first case, whereas in the 
second we evade the problem by insisting that on physi- 
cal grounds one must average the cross-section for emission 
of soft photons (Brehmsstrahlung) over some experimental 
resolution for the detection of the photon energy™ 

The crucial observation is that while the Green's functions 
may have IR divergences the cross-sections smeared over ex- 
perimental resolutions are always finite. At the diagram- 
matic level what happens is that there are cancellations be- 
tween loop integrals (or Feynman parameter integrals) and 
phase space integrals, so our intention is to extend our pro- 
gram to treat loop and phases space integrals consistently. 
What we must do is to ensure that all divergences cancel lo- 
cally in the integrand, just as we do for Green's functions by 
our use of the R operation. To this end we intend to gener- 
ate the integrands for cross sections by representing them as 
cut bubble diagrams, where the cut propagators correspond 
to on-shell particles convoluted with the experimental reso- 
lution functions. 

9. COMPUTATIONAL ISSUES 

Our prototype program is written in Maple and generates 
C code: an example of its use is illustrated in Figures [1] 
and [2] This is unusual in that most symbolic computations 
of Feynman diagrams make use of domain-specific languages 
such as FORM 42 , 47 , 48 , which are relatively unsophisticated 



8 At least "if we keep our feet on the ground and ignore 
gravity" (E. Mottola). 

9 In fact probably only the photon, as neutrinos seem to have 
masses, albeit very small ones. Gluons, the gauge bosons for 
QCD, are also massless; but as QCD is a strongly interact- 
ing confined theory they do not exist as physical "on-shell" 
particles. 

10 In mathematical terms we say that the Green's functions 
in Minkowski space are (tempered) distributions rather than 
functions. 



diag2 := diagram ( 

seq(vertex(cat ( ' V , i) , 'ThreeScalarVertex' ) , i = 4.. 7), 

externalline (E4, V4,p3 , ' Scalar ' ) , externalline (E5 , V6 , -p3, 'Scalar' ) , 

line (14, V4.V5, 'Scalar') , line(I5,V4,V7, 'Scalar') , line(I6,V5,V7, 'Scalar') , 

line (17, V5.V6, 'Scalar') , line(I8,V6,V7, 'Scalar'))! 

g2 := analyzeDiagram(diag2) $ 

analyzeDiagram: "2 loop diagram" 

rg2 := R(g2); 

R: "Applying R to IPI graph Graph(EXT(E4,E5) , INT(I4, 15, 16, 17 , 18) ) " 
R: "Applying R to IPI graph Graph(EXT(E4, 17 , 18) , INT (14, 15, 16) ) " 

T: "Taylor subtracting 1 loop graph Graph(EXT(E4, 17 , 18) , INT(I4, 15 , 16) ) of degree 0" 
R: "Applying R to IPI graph Graph(EXT(E4,E5 , 16) , INT (14, 15, 17, 18) ) " 

T: "Taylor subtracting 1 loop graph Graph(EXT(E4,E5 , 16) , INT (14, 15 , 17 , 18) ) of degree -2" 
R: "Applying R to IPI graph Graph(EXT(E5, 14, 15) , INT(I6, 17, 18) ) " 

T: "Taylor subtracting 1 loop graph Graph(EXT(E5, 14, 15) , INT(I6, 17 , 18) ) of degree 0" 
T: "Taylor subtracting 2 loop graph Graph(EXT(E4,E5) , INT (14, 15, 16 , 17, 18) ) of degree 2" 
T: "Taylor subtracting 2 loop graph -Graph (EXT (E4.E5) , INT (14, 15 , 16 , 17 , 18) ) of degree 2" 
T: "Taylor subtracting 2 loop graph -Graph (EXT (E4.E5) , INT (14, 15 , 16, 17, 18) ) of degree 2" 



Figure 1: Example of the use of our Maple program for the two-loop two-point function 4 \ in 4> A theory. 



but can execute their limited repertoire of operations on very 
large expressions with great efficiency. 

9.1 Code Generation 

We could try to evaluate the final numerical integrals 
within Maple, but other than for debugging purposes this is 
too inefficient even using the NAG integration routines pack- 
aged within Maple. We therefore need to generate efficient 
numerical code to evaluate the integrand for use within our 
special-purpose Monte Carlo integration programs. It is not 
too important which language is used, our example is in C 
but it would be trivial to generate Fortran if that was pre- 
ferred. 

9.2 Memory Management and Laziness 

The computational model that FORM uses is to stream a 
sum of many terms from disk to disk, either applying pat- 
tern-based transformations or sorting and collecting terms 
as the data passes through the processor. This model is jus- 
tified by the huge size of the expressions that are generated 
in realistic Feynman diagram computations, unlike the toy 
example of Figure [T] While our Maple program works well 
on moderate size problems it dies a horrible death when the 
entire problem no longer fits in memory. 

We suggest that the huge intermediate expressions gener- 
ated are not really an intrinsic property of the problem, but 
are an artefact of the way we think about and implement 
the programs. To be specific, our program first applies the 
R operation generating a sum of Taylor subtraction terms 
as it recurses through the Henge decomposition; it then runs 
through all these terms converting them into Feynman pa- 
rameter integrals; it simplifies the tensor structures in the 
numerators; and it then collects all the terms with similar 
quadratic forms in their denominator before finally gener- 
ating C code. This multi-pass approach is not necessary, 
it is just convenient for thinking about and debugging the 
program, as it separates the computation into logically in- 
dependent steps. 



We could avoid generating huge intermediate expressions 
by generating the terms lazily, and this could be done declar- 
atively without altering the logic of the program if Maple im- 
plemented a lazy map primitive, or a yield statement for use 
within loops. We would then generate each term and imme- 
diately apply subsequent operations on it until we reach the 
"collection" phase, in which terms could be collected into 
elements of a hashed array. 

Eventually this approach will still run out of memory if 
this array has to live in (virtual) memory, but this could 
be circumvented by writing the array elements to disk when 
necessary. 

10. CONCLUSIONS 

Current perturbative calculations in QFT are unthinkable 
without the use of computer algebra. Whether this requires 
domain-specific languages such as FORM or more careful mem- 
ory choreography in general purpose systems such as Maple 
remains to be seen. All of the pieces of the calculation, in- 
cluding renormalization, can be fully automated except for 
evaluation of the Feynman parameter integrals where the 
lack of closed-form solutions forces us into the arms of nu- 
merical integration. 
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